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1.  INTRODUCTION  AND  SUMMARY 


Recent  progress  in  the  development  of  high  power 
lasers  has  evolved  lasers  capable  of  causing  thermal  dam- 
age to  material  surfaces.  In  order  to  adequately  evaluate 
the  effects  of  material  properties  tnd  laser  beam  charac- 
teristics on  material  damage,  a reliable  analytical  model 
is  needed.  After  reviewing  many  of  the  available  analy- 
tical models,  it  was  found  that  most  were  limited  to 
one-dimensional  heat  conduction,  constant  thermal  proper- 
ties, constant  irradiation,  and  no  radiation  relief  or  aero- 
dynamic cooling.  In  many  flight  applications  these  effects 
are  significant.  Therefore,  extensive  modifications  have 
been  made  to  the  APL  finite  difference  heat  transfer  pro- 
gram (Ref.  1)  to  analyze  laser  heating  of  flight  vehicles. 

This  new  program  is  known  as  the  Continuous  Wave  Laser 
Damage  Computer  Program  (CLAD).  It  computes  mass 
loss  and  temperature  histories  of  laser  heated  bodies  ac- 
counting for  three-dimensional  conduction,  temperature  de- 
pendent thermal  properties,  vaporization,  melting,  chemi- 
cal reactions,  aerodynamic  heating,  radiation  relief,  and 
material  removal. 

The  APL/BBE  Continuous  Wave  Laser  Damage  Com- 
puter Program  (CLAD)  uses  finite  difference  techniques  to 
predict  temperature  histories  in  Icser  heated  materials. 

The  program  is  written  in  PL/I  computer  language,  and  in- 
cludes: (a)  three-dimensional  conduction,  (b)  thermal 
properties  as  functions  of  temperature  (thermal  conduc- 
tivity, thermal  capacitance,  absorptivity,  and  emissivitj  \ 
(c)  radiation  relief,  (d)  latent  heat  of  fusion,  and  (e)  laser 
irradiation  of  each  surface  element  specified  as  a function 


Ref.  1.  D.  W.  Fox,  H.  Shaw,  and  J.  Jellinek,  "Nu 
merical  Approximations  in  Heat  Transfer  Problems  and 
Usage  of  IBM  7090  Computer  for  Solutions,  " ATL/JHU 
CF-2954,  17  May  1962. 
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of  time,  with  on  option  to  specify  irradiation  by  a moving 
Qausaian  beam.  CLAD  haa  recently  been  expanded  to  in- 
clude: (a)  aerodynamic  heating  and  cooling,  (b)  material 
vaporiaation , fc)  chemical  reactiona  at  the  surface, 

<c)  material  removal  of  vaporised  and  ablated  elements, 
and  (e)  cylindrical  elements. 

This  report  has  been  written  as  a user’s  manual  for 
the  CLAD  program  and  includes  an  overview  of  the  program, 
a discussion  of  assumptions  and  limitations,  a main  pro- 
gram description  and  a sample  problem. 
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2.  OVERVIEW  OF  CLAD  CALCULATIONS 


The  following  discussion  presents  a brief  overview 
of  calculations  performed  in  CLAD,  and  is  followed  by  de- 
tailed discussions  of:  (a)  the  surface  energy  balance, 

<b>  time  step  calculations,  <c)  corrections  for  blowing, 

<d)  an  option  for  mechanical  erosion,  and  (e)  assumptions 
and  limitations. 

The  CLAD  computer  program  is  a finite  difference 
program  written  for  analysis  of  laser  heated  plates  and 
cylinders.  The  program  includes  radiative  and  convective 
heating  with  heat  relief  via  radiation  and  conduction.  In 
addition,  the  program  can  account  for  in-depth  melt.’ng  and 
material  removal  at  the  heated  surface.  The  material  re- 
moval process  can  be  via  vaporisation  or  chemical  reac- 
tions at  the  surface.  The  vaporization  process  occurs  at  a 
specific  temperature  and  heat  of  vaporization,  while  the 
chemical  reactions  occur  over  a range  of  surface  tempera- 
tures and  either  absorb  or  release  heat  at  a rate  dependent 
on  the  surface  temperature  and  the  chemical  nature  of  the 
ablating  materiel. 

A pictorial  view  of  some  of  the  parameters  included 
in  CLAD  is  presented  in  Fig.  1,  The  program  commutes  a 
surface  temperature  (Tw)  of  each  zero  volume  surface  node, 
by  iteratively  solving  an  energy  balance  at  the  surface.  For 
ablating  surfaces  the  surface  energy  balance  includes  w 
chemical  heatir.g  term  which  is  obtained  from  the  Equilib- 
rium Surface  Thermochemistry  Program  (EST)  (Ref.  2),  as 
a function  of  surface  temperature  and  local  pressure.  Fol- 
lowing computation  of  the  surface  temperature,  CLAD  com- 
putes the  ablation  mass  loss  rate  as  a function  of  the  surface 
temperature  and  local  pressure,  Capacitances  of  all  finite 
volume  elements  (2  and  3)  are  then  computed  as  functions  of 
temperature,  and  conduction  heating  rates  are  computed 

Ref.  2.  User's  Manual,  Aerotherm  Equilibrium 
Surface  Thermochemistry  Program,  " Version  3,  Aero- 
therm  Corporation,  Mountain  View,  California,  ! eport 
UM-70-13,  May  1970, 
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between  adjacent  nodes.  A critical  stability  time  strp  ia 
then  used  to  compute  node  temperatures  at  the  next  time 
step.  Corresponding  surface  temperatures  are  then  com- 
puted for  the  sero  volume  surface  elements*  and  the  above 
calculations  are  repeated  for  the  next  time  step.  This  se- 
quence continues  until  the  end  time  is  reached. 


SURFACE  ENERGY  BALANCE 

The  chemical  reaction  model  used  in  CLAD  is  the 
same  at*  used  in  Ref.  3.  In  thid  model  (Fig.  1)  the  surface 
energy  balance  per  unit  area  is; 

q.  y q , + q q . - q . 4 - q *0.  (1) 

in  xhem  x:onv  -cond  ^rad  out  ^sens 

Looking  at  the  individual  terms,  the  incoming  heat 
rate  (qin)  is  made  up  of  laser  heating  and  radiative  heating 
as  follows; 


B 

B 

8 

I 

ii 

1 


'in 


|(qlasl  (view)  . 


t 


4rai  in 


ALPHA  , 
w 


where: 


(2) 


^las 

VIEW 

**rad  in 


is  the  incident  laser  heat  flur  per  unit 
area, 

is  a view’  factor  to  account  for  angular 
beam  impingement, 

is  an  incident  radiation  heat  flux  per 
unit  area,  and 


ALPHA  is  the  surface  absorptivity  evaluated  at 

W T 

1 W* 


Ref.  3,  "User's  Manual  Aerotherm  Charring  Mate- 
rial  Thermal  Response  and  Ablation  Program,"  Version  3, 
Aerotherm  Corporation,  Mountain  View,  California,  Re- 
port RM-70-14,  April  1070. 
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Like  all  uuch  models*  the  model  used  to  obtain 
qChem  is  *n  approximation.  A discussion  of  its  validity  is 
beyond  the  scope  of  this  report  but  has  been  made  in  Refs. 

4 and  5.  In  thie  model  (Fig.  1)  free  stream  gases  diffuse 
through  the  boundary  layer  to  the  surface,  where  they  react 
with  the  ablating  material  and  either  release  or  absorb  heat. 
The  resulting  gases  are  then  removed  from  the  surface  via 
diffusion  and  convective  mass  transfer  into  the  boundary 
layer.  Appendix  A presents  a derivation  of  the  chemical 
heating  rate  (qchem)  in  terms  of  diffusion  and  mass  trans- 
fer rates. 

However,  the  chemical  reaction  model  contains  sev- 
eral important  assumptions,  which  are: 

1.  All  chemical  reactions  take  place  at  the  wall 
and  no  reactions  occur  in  the  boundary  layer. 

2.  No  solid  reactants  are  produced  at  the  wall  and 
no  ablative  material  leaves  the  surface  due  to 
mechanical  erosion;  For  special  cases  where 
ATJ  graphite  is  used,  mechanical  erosion  can 
be  included  (see  Mechanical  Erosion  Option 
later  in  this  Section). 

3.  There  is  no  pyrolysis  within  the  ablative  mate- 
rial; all  pyrolysis  is  confined  to  the  surface. 


Ref.  4.  C.  B.  Moyer  and  R,  A.  Rindall,  "Finite  Dif- 
ference Solution  for  the  In-Depth  Response  of  Charring 
Materials  Considering  Surface  Chemical  and  Energy  Bal- 
ances, " Aerotherm  Corporation  Final  Report  66-7,  Part  II, 
14  March  1967  (also  NASA  CR-106I,  June  1968). 

Ref,  5,  E,  P,  Bartlett  and  R.  D.  Grose,  "An  Eval- 
uation of  a Transfer  Coefficient  Approach  for  Diffusion  Coef- 
iicients,  " Aerotherm  Corporation,  Mountain  View,  Calif- 
ornia, Report  69-50,  30  June  1969. 
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4.  Gases  diffuse  through  the  boundary  layer  accord' 
ing  to  their  mass  concentrating  gradients  and 
these  gradients  are  assumed  linear  through  the 
boundary  layer. 

5.  There  is  no  thermal  diffusion  of  gases  in  the 
boundary  layer. 

The  chemical  heating  term  (qcftem)  °*  Etl*  is 
made  up  of  two  parts.  The  first  term  (q^)  is  the  rate  of 
energy  release  due  to  diffusing  gases  reacting  at  the  sur- 
face. In  this  process,  each  free  stream  specie  i diffuses 
through  the  boundary  layer  to  the  surface  where  it  reacts 
with  the  surface  and  releases  heat.  The  product  specie  i 
then  leaves  the  surface  at  a total  enthalpy  (AHit  w)  different 
than  it  had  on  arrival.  This  heating  rate  is  expressed  as: 


= p U 
e e 


N 

<c  >h  E 

mbi=l 


(Z.  - Z.  ) AH. 

x,  e x,  w i,  w 


(3) 


where: 


P U 
e e 


N 

(C  ).  E (Z.  - Z.  ) 

m b^=1  i,e  x,  w 


is  the  diffusion 


mass  flow  rate  of  specie  i arriving  at  a 
unit  surface  area. 


(C  is  the  mass  transfer 

Stanton  number  (md/reUeAZ)|j,  corrected 
for  blowing  as  shown  later  in  this  Section 
(Blowing  Rate  Correction),  with  AZ  = the 
diffusion  driving  potential  which  produces 
the  diffusion  mass  flow  rate  (m^)  into  the 
surface. 


Z,  and  Z.  are  the  diffusion 

x e x,  w 

driving  potentials  for  the  ith  specie  diffus 
ing  through  the  boundary  layer  to  the 
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surface  (Zjt  e is  the  potential  at  the  bound- 
ary layer  edge  and  Zjt  w is  the  potential  at 
the  wall),  and 

is  the  enthalpy  change  of  each  diffusing 
' specie  due  to  its  reaction  with  the  sur- 

face material  (AHjtW  = Hj>w  - Hrg  W). 

The  second  part  of  the  qchem  term  represents  the 
rate  of  energy  released  (qa)  by  the  ablating  material  as  it 
changes  from  a solid  ablative  to  reacted  gases: 


q =p  U (C  ).  $ 

"o  O A K 


m 


(4) 


where 


p U (C  ).  j9  is  the  ablation 

e e mb 

rate  (ma) 

0 is  the  ablation  to 

diffusion  mass  flux  ratio  defined  as 
3 - nia/peU0 

H is  the  total  enthalpy 

9 (above  Tref)  of  the  solid  ablative  material 

evaluated  at  the  wall  temperature  (Tw),  and 

H is  the  total  enthalpy 

6*  (above  Tref)  of  the  gaseous  reactants  eval- 
uated at  Tw  (HrgfW  is  the  sum  of  the  chemi- 
cal and  sensible  enthalpy  of  the  gases  and 
is  averaged  according  to  the  mass  fraction 
of  each  gas  present  at  the  surface). 

Hence,  the  total  chemical  heating  term  is: 


^chem  + 


= p U 

e e 


<C  ), 
m o 


N 

E 

i- 1 


(Zi.e 


- Z.  )AH.  + /S(H 
i,  w i,  w i 


a,  w 


H ) 
rg,  w 


. (5) 
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In  Eq.  (1)  the  convective  heat  transfer  term  is: 

q ^ p U (C-J.  (H  - H ) . 

Mconv  e e H b r e,  w 


where: 


PeUe  ^CH)b 


is  the  enthalpy  based  heat  transfer 
coefficient  corrected  for  blowing 
as  will  be  discussed  later. 


is  the  heat  transfer  Stanton  number 
corrected  for  blowing, 

H is  the  total  recovery  enthalpy 

1 (above  Tref^’  and 

H is  the  total  enthalpy  (above  Tref) 

e* w of  the  gases  at  the  boundary  layer 

edge  and  evaluated  at  the  wall  tem- 
perature. 

The  convective  heat  flux  is  more  usually  written  in 
the  form: 

q = p U (C„).  (h  - h ) , (7) 

^conv  e e H b r w 

where : 


h and  h 
r w 


are  the  sensible  recovery  and  wall 
enthalpies,  respectively. 


The  derivation  of  Eq.  (6)  from  Eq.  (7)  is  given  in  Refs.  4 
and  6. 


Ref.  6.  R,  M.  Kendall,  R.  A.  Rindall,  and  E.  P. 

Bartlett,  "A  Multicomponent  Boundary  Layer  Chemically 
Coupled  to  an  Ablating  Surface, " AIAA  Journal,  Vol.  5, 

No.  6,  June  1967,  pp.  1063-1071. 
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For  unity  Lewis  number,  qchem  and  ^conv  can 
combined  and  simplified.  When  the  Lewis  number  equals  1, 
the  mass  and  thermal  boundary  layers  are  coincident  in 
s.pace  and  Z|(  e is  evaluated  at  the  thermal  boundary  layer 
edge.  Hence: 

(CH>b  = <Cm'b  • 


<Cm>b 


' Zi,  w*  iHl.w  5 <CH*b 


■ 

J ■< 

H - H 

• 1 

e,  w rg,  w 

l\ 

t ! ■$ 

For  the  unity  Lewis  number: 


q + q . = p U (CTT), 

Mconv  Mchem  e e H b 


(H  - H ) + )3  (H  - H ) 
r rg,  w a,  w rg,  w 


where  the  heat  transfer  coefficient  peUe  (CH)  and  recovery 
enthalpy  Hr  are  time  dependent  variables  input  by  the  user. 
The  heat  transfer  coefficient  is  automatically  adjusted  for 
blowing  as  is  discussed  later  in  this  Section.  The  quantity 
Ha  w is  obtained  from  a utility  program  which  computes 
Hft  w as  a function  of  wall  temperature  as  follows; 


w (pc  ) 

-5  -7^-  dT  • 


The  remaining  terms  in  Eq.  (8)  are  obtained  from 
the  Equilibrium  Surface  Thermochemistry  Program  (EST) 
(Ref.  2).  For  unity  Lewis  number  EST  supplies  Tw  and 
Hr„  w as  tabular  functions  of  local  pressure  (PL)  and  the 
ablation  to  diffusion  ratio  (/3).  As  part  of  the  CLAD  pro- 
gram, a utility  subroutine  accepts  EST  output  directly  and 
converts  it  to  tables  of  HrgtW,  j8  (Ha#  w - Hrg,  w),  and  j3  as 
functions  of  Pl  and  T. 

For  Lewis  numbers  not  equal  to  1,  qcfcem  anc*  qconv 
are  as  presented  in  Eqs.  (5)  and  (6),  and  the  EST  program 
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supplies  Tw#  HrgjW,  HGj  w,  and  Z (Zj>e  - Z^w)  AH^W 

as  tabular  functions  of  PL  and  0.  The  utility  subroutine 
accepts  this  data  directly  and  converts  it  to  tables  of 

V »•  «e.W  t*  <Zi.  e - zi.  w>  iHi.  w H.g,  „>• 

and  0 as  functions  of  PL  and  Tw. 


The  remaining  terms  in  Eq.  (1)  are: 


^cond 

w-2 
L _ 

<T2 

- V- 

w-2 

^rad  out 

= Fafw 

(T4 

w 

- T<p>.  and 

q 

- m (n 

h 0)  . 

^sens 

a a,  w 

a,  2 

where: 

K is  the  material  thermal  conductivity  eval- 

w uated  at  (Tw  + T2)/2, 

L.  _ is  the  conduction  path  length  from  the  wall 
W to  node  2, 

F is  a view  factor, 

a is  the  Stefan- Boltzman  constant. 


c is  the  surface  emissivity,  and 

W 

q is  the  sensible  energy  required  to  raise 

the  ablating  material  from  temperature 
T2  to  surface  temperature  Tw. 

In  Ref.  3 this  is  referred  to  as  the  "velocity  term. " 
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TIME  STEP  CALCULATIONS 


At  the  beginning  of  time  step  a , Eq.  (1)  is  solved 
iteratively  for  surface  temperatures  at  time  ra.  For  ablat- 
ing surfaces  the  ablation  to  diffusion  ratio  (fi)  is  found  from 
tables  as  a function  of  Tw  and  PL.  fi  is  then  used  to  find 
the  ablation  mass  loss  rate  per  unit  area  from: 


’ * (<e*Ue  <Cm*b] 


where: 


P U <C  ).  = P U (CJ.  (C  / C__)  . 
ee  mb  ee  Hb  m H 


and  PeUe  Cjj  is  input  as  a function  of  time  and  corrected  for 
blowing  as  will  be  shown  later.  Cm/Cpj  is  also  an  input 
value. 


For  surfaces  that  are  vaporizing  instead  of  chemi- 
cally ablating,  the  surface  energy  balance  becomes; 

q + q ■+•  m Q +q  - q , >•  q =0,  (10) 

Min  Mconv  a ^vap  Mcond  Mrad  out  Msens 

where  Qvap  Is  the  heat  of  vaporization  per  pound  of  mate- 
rial. By  setting  the  surface  temperature  equal  to  the  vapor 
temperature,  Eq.  ( 10)  is  used  to  compute  ma. 

Following  computation  of  ma»  the  CLAD  program 
computes  heat  conduction  between  elements  and  elemental 
capacitances  These  values  are  then  used  in  the  Dusinberre 
stability  equation  (Ref.  7 ) to  find  the  stability  time  step  which 
equals  90%  of  the  smallest  value  of: 


<P°  v>i 


Ref.  7.  G.  M.  Dusinberre,  Heat  Transfer  Calcula- 
tions by  Finite  Differences.  International  Textbook  Co. , 
Scranton,  1961. 
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where  i is  all  elements  having  a capacitance  associated 
with  them,  and  the  sum  is  taken  over  all  elements  j con- 
nected to  i.  The  temperature  of  interior  elements  2,  3, 
etc.  are  then  computed  at  time  If  any  element  has  gone 
above  its  melt  temperature,  then  subroutine  QFUS3N  com- 
putes the  melt  depth  in  each  element.  If  an  element  has  not 
completely  melted,  then  it  is  set  to  the  melt  temperature. 

If  the  element  is  resolidifying  then  QFUSON  accounts  for 
the  heat  of  fusion  and  reduces  the  melt  depth  accordingly. 

The  laser  surface  heating  rate  (Qi*8)  Is  next  com- 
puted at  time  r^.  CLAD  then  uses  the  mass  loss  rate  com- 
puted at  t4  to  determine  the  volume  of  element  2 removed 
during  the  time  interval  - ra.  As  the  volume  of  element 
2 decreases,  the  stability  time  step  approaches  zero.  There- 
fore, after  a specified  portion  of  element  2 (usually  50%)  has 
ablated  away,  the  program  lumps  the  remaining  portion  with 
element  3.  Element  2 then  becomes  a zero  volume  surface 
node  and  element  1 is  removed  from  the  model  (as  shown  at 
Tc  in  Pig.  2).  In  this  manner  elements  are  removed  with- 
out seriously  reducing  the  stability  time  step. 

As  each  element  shrinks,  conduction  paths  to  adja- 
cent elements  are  also  changed.  Conduction  paths  perpen- 
dicular to  the  surface  are  adjusted  according  to  the  distance 
between  elements,  assuming  the  node  always  remains  in  the 
center  of  each  element.  It  is  important  to  remember,  that 
during  mass  removal,  the  location  of  element  2 is  continu- 
ally moving  toward  the  backface. 

Thermal  resistance  between  longitudinally  adjacent 
elements  are  computed  as  the  sum  of  the  resistances  of 
each  element.  For  example  the  thermal  resistance  between 
elements  2 and  5 in  Fig.  2 at  is; 
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between  elements  3 end  6 Is  calculated  In  the  seme  manner 
as  between  2 and  5 in  Eq.  (11).  Notice  that  even  though 
elements  3 and  5 have  a common  Interface  they  are  not  di- 
rectly connected;  however,  they  are  indirectly  connected 
through  element  6. 

After  computing  the  thermal  resistance  between 
ablating  elements  the  program  calculates  surface  tempera- 
tures at  time  TV  using  Eq.  (1),  as  done  previously  at  time 
r . The  procedure  is  continued  computing  ma,  T2.  Tg, 
etc.  at  each  time  step  until  the  end  time  is  reached. 


BLOWING  RATE  CORRECTION 


The  CLAD  program  automatically  reduces  the  heat 
and  mass  transfer  Stanton  numbers  (C^  and  Cm)  to  account 
for  the  familiar  blowing  effect  according  to  the  equation: 


^CH^b  = CH  Ot 


e - 1 


where: 


2Xm 


*eUe  CH 


(12) 


(13) 


and  X is  an  empirical  constant  known  as  the  "blowing  parame- 
ter. " Reference  3 recommends  that  X = 0. 5 be  used  for  lami- 
nar flow  and  that  X = 0.  4 correlates  with  turbulent  data  some- 
what better. 


The  CLAD  program  computes  (C^^  from  an  analogy 
bewteen  heat  and  mass  transfer,  where  (Cm)jj  = (Cjj)^  Cm/CH 
and  the  cor«*tant  ratio  Cm/Cfj  is  input  to  the  program. 
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For  graphlti,  the  user  has  the  option  of  using  the 
Puts  and  Bartlett  correlation  (Kefs.  8 and  8)  instead  of 
Eq.  (18).  The  Puts  and  Bartlett  equations  are  as  follows: 


<CU>b  * CH  [ * 0,°lf94  •*  ♦ O.OSStS  I*  < 0,0)115  |*j  , (14) 

and 


«Vb 


H 


- 1 


(15) 


where: 


;r  ) 

' "m  b 


(16) 


and 


X « (0. 012  + 0. 018  8 + 0. 0814  fl2)  * (0. 762  - 0. 038  8 ).  (17) 

o o o 


The  user  . .'tiat^s  the  Puts  and  Bartlett  correlation  for 
graphites  by  changing  CMjCH  to  a negative  value  in  the 
CALL  ivBAROS  statement  (see  Section  3).  This  also  ini'* 
tiates  the  mechanical  erosion  option  discussed  below. 


MECHANICAL  EROSION  OPTION 

In  the  above  discussion  of  CLAD,  the  ablation  mass 
loss  (ma)  is  computed  assuming  no  mechanical  erosion.  At 

J^TT  K.  E.  Puts  and  E.  P.  Bartlett.  "Heal  Trans  - 
fer  and  Ablation-Rate  Correlations  for  Reentry  Heat  Shield 
and  Nose  Tip  Applications,"  AIAA  Preprint  Paper  No,  72-91, 
January  1974. 

Ref,  9,  L.  L.  Perini,  "Heat  and  Maos  Transfer  Cor- 
relation Equations  for  Subliming  Graphite  in  High  Speed 
Flow,  ' APL/JHU  ANSP-M-11,  August  J 974. 
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high  l*nv«i*iures  and  in  a high  shear  environment,  mechani- 
cal erosion  can  became  significant.  Although  CLAD,  in  gen- 
eral. cainot  handle  mechanical  erosion,  it  does  have  a spec- 
cial  calculation  for  ATJ  graphite,  where  the  total  mass  loss 
rate  is  the  sum  of  the  chemical  ablation  mass  loss  rate  (ma) 
and  the  mechanical  mass  loss  rate  (mmec^>.  In  Ref.  10  the 
mmech  *8  derived  from  the  difference  between  experimentally 
measured  mass  loss  data  of  Lundell  and  Dickey  (Ref.  1 1 ) and 
analytical  values  computed  neglecting  mechanical  erosion.  A 
curve  through  the  bulk  of  the  points  gives: 

"Wh  ' m»  [ .1*76  -Ttw/1.8)  ' 155]  • <18> 

which  is  applicable  over  the  temperature  range  4500  * T * 
7245°R  and  local  stagnation  pressure  range  0,  3 * P^,  * 4.A 
atm.  Since  no  experimental  data  are  available  above  7245°F. 
mmech  is  set  ®duai  to  J*89  ma  at  values  above  7245°R.  It 
should  be  emphasised  that  the  above  calculation  is  only  valid 
for  ATJ  graphite.  The  user  may  initiate  the  mechanical  ero- 
sion option  by  changing  CMjCH  to  a negative  number  in  the 
CALL  ABAROS  statement.  ~ 


CLAD  - ASSUMPTIONS  AND  LIMITATIONS 

There  are  several  assumptions  and  limitations  in 
CLAD  of  which  the  user  should  be  aware.  These  are  dis- 
cussed in  the  following  paragraphs: 

1.  The  CLAD  program  is  limited  to  only  one  ablat- 
ing material  in  any  given  analysis.  This  restric- 
tion can  be  eliminated  in  the  future  if  the  need 
arises,  but  will  require  some  reprogramming. 

Ref.  To!  L.  L.  Perini,  "Review  of  Graphite  Ablation 
Theory  and  Experimental  Data,"  APL/JHU  ANSP-M-1, 
December  1971. 

Ref,  11.  J.  H.  Lundell  and  R,  R.  Dickey,  "Graphite 
Ablation  at  high  Temperatures, " AIAA  Preprint  Paper  No. 
71-418,  April  1971. 
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Although  CLAD  ia  limited  to  one  ablating  mate- 
rial there  may  bo  several  substrates  of  differ- 
ent materiala  aa  long  aa  theca  aubatratea  don't 
ablate.  In  addition  CLAD  ia  restricted  to  models 
which  either  chemically  ablate  or  vaporise  but 
cannot  handle  modelo  which  do  both.  For  one- 
dimensional  models,  both  these  restrictions  can 
be  handled  by  running  the  program  twice.  The 
first  run  was  continuous  until  the  first  surface 
layer  is  completely  ablated,  then  the  second  run 
is  made  using  the  second  surface  material  ini- 
tialised to  the  temperatures  at  the  end  of  the 
first  run. 

The  user  should  realise  that  although  CLAD 
ia  restricted  to  a single  ablative  per  run.  it  is 
not  restricted  to  a single  vaporising  material 
per  run.  It  can  have  several  material  layers 
each  vaporising  when  it  becomes  the  surface 
layer. 

2.  As  elements  ablate  or  vaporise  the  surface  tem- 
perature fluctuates  in  a "sawtooth"  manner 
caused  by  node  lumping.  This  sawtoothing  is 
usually  not  too  severe  and  may  be  reduced  by 
using  smaller  elements  (Ref.  J2). 

3.  If  during  a single  time  step,  a vaporising  sur- 
face has  not  yet  melted,  then  its  heat  of  fusion 
is  ignored  for  that  time  step.  This  is  usually 
rot  very  serious  since  the  heat  of  vaporisation 
is  usually  much  larger  than  the  heat  of  fusion. 

4.  As  surface  elements  ablate  the  node  adjacent  to 
the  surface  is  always  located  in  the  center  of  the 
element.  Hence  as  the  element  ablates  the  node 
is  continually  moving  towards  the  backface.  This 
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can  sometimes  be  confuting  whan  plotting  a tem- 
paratura  history  of  this  node,  tinea  tha  ‘ lament 
location  ia  not  constant. 

5.  No  attempt  haa  bean  made  to  haat  tha  aida  of 
thoaa  elements  whoaa  neighboring  elements  have 
ablated  or  vaporised. 

6.  Masa  removal,  chemical  reactions,  and  vapori- 
sation occur  only  at  tha  aurfaca. 

7.  The  number  of  different  materials  is  presently 
limited  to  nine.  If  additional  materials  are  re- 
quired. CLAD  could  be  reprogrammed  to  handle 
them. 

Although  the  CLAD  analysis  contains  several  assump- 
tions and  limitations,  it  is  still  quite  general  and  can  be  uti- 
lised for  a wide  variety  of  applications.  The  programming 
details  necessary  for  implementing  CLAD  are  discussed  in 
the  following  section. 


3.  MAIN  PROGRAM  DESCRIPTION 
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The  CLAD  program  consists  of  a aeries  of  subrou- 
tines called  by  a main  program.  A typical  main  program 
is  as  follows: 


NAME:  PROCEDURE  OPTIONS  (MAIN) 

DECLARE  STATEMENTS 
CALL  STORE 

A series  of  initialisation  statements 
which  might  include:  GET  COPY 
DATA,  GET  LIST.  CALL  READRK, 
CALL  READTM,  CALL  READAB. 
and  CALL  RFnCP3 

CALL  THREED* 

CALL  SET 

11:  IMAX  = DECIDE  (TIMS,  IOX,  TIMEI)' 

CALL  GAUSBMV 
CALL  QFUSON' 

CALL  ABAR03 
CALL  CAPCN3 
CALL  WRITE 
CALL  STEP 


♦These  statements  are  optional. 


PrinJia!  MU  link 
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END  NAME 

The  order  in  which  the  above  subroutines  are  called 
is  very  important  and  should  be  exactly  as  shown  for  the 
program  to  function  correctly. 

A complete  description  of  each  section  of  the  main 
program  is  presented  below.  In  reading  this  material,  the 
reader  may  find  it  helpful  to  refer  to  Appendix  B,  which 
presents  the  main  program  used  to  analyze  a sample  prob- 
lem discussed  later  in  the  text. 

1.  NAME:  PROCEDURE  OPTIONS  (MAIN)  de- 

fines procedure  NAME  as  the  main  program. 

2.  DECLARE  STATEMENTS  define  all  variables 
and  subroutines  used  in  the  main  program.  All  subroutines 
are  declared  as  entry  points  using  DCL  ENTRY  statements 
exactly  as  they  appear  in  the  sample  problem  in  Appendix  B. 

In  addition  all  variable  names  used  in  the  program  are  de- 
clared with  the  following  attributes:  FIXED  or  FLOAT; 

BINARY  or  DECIMAL;  AUTOMATIC,  STATIC,  CONTROLLED, 
or  EXTERNAL;  and  double  or  single  precision.  The  following 
variable  names  must  be  declared  EXTERNAL  and  must  always 
appear  in  the  declare  statements  QVAPOR(*),  TVAP(*)f  RO(*), 
and  AIRPRNT.  Additional  variables  declared  EXTERNAL  for 
GAUSBM  or  QFUSON  are:  IMAX,  QMELT(*),  TIMS,  TMELT(*), 
and  NDX  (20,20).  In  addition  variable  names  QVAPOR,  TVAP, 
QMELT,  TMELT,  and  RO  are  specified  as  CONTROLLED  and 
must  be  assigned  array  lengths  in  an  ALLOCATE  statement. 

The  array  length  is  equal  to  the  number  of  vaporizing  mate- 
rials. These  variable  names  are  defined  as  follows: 

QVAPOR  is  an  array  defining  the  latent 

heat  of  vaporization  for  each 
material  assigned  thermal 
properties  in  READRK,  where 
the  subscript  of  QVAPOR  is 
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the  index  defining  the  mate- 
rial  (FLOAT  BINARY) 

(Btu/lb), 

TVAP  is  an  array  of  corresponding 

vapor  temperatures  (FLOAT 
BINARY)  (°R);  TVAP(»)  is 
zero  if  material  I ablates 
with  a chemical  reaction  in- 
stead of  vaporizing*  (In  addi- 
tion, for  cases  with  vaporiza- 
tion, Hpg  is  set  to  air  enthalpy 
values  in  the  input  to  READAB. 
Hence  CLAD  is  restricted  to 
models  which  either  vaporize 
or  chemically  ablate,  and  can- 
not handle  models  with  both), 

QMELT  is  an  array  defining  the  latent 

heat  of  fusion  for  each  mate- 
rial assigned  thermal  proper- 
ties in  READRK,  where  the 
subscript  of  QMELT  is  the  in  - 
dex  defining  the  material 
(FLOAT  BINARY)  (Btu/lb), 

TMELT  is  an  array  of  corresponding 

melt  temperatures  (FLOAT 
BINARY)  (°R),  and 

RO  is  an  array  of  corresponding 

densities  (FLOAT  BINARY) 

(lb /ft3). 

The  attributes  assigned  to  other  variables  in  the  declare 
statements  are  discussed  below  as  they  are  encountered  in 
the  program. 


3.  CALL  STORE  (LASCAP,  LASCAP,  3)  specifies 
all  ailocatable  variables  as  arrays  of  length  LASCAP,  and 
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initializes  them  to  0. 0.  LASCAP  is  a fixed  binary  number 
equal  to  the  largest  capacitor  number  used  in  the  program. 
The  ”3"  informs  STORE  to  allocate  variables  used  in  rou- 
tines THREED,  QFUSON,  ABAR03.  and  CAPCN3. 

4.  If  variable  names  are  not  assigned  initial 
values  in  the  above  DECLARE  statements,  they  are  given 
values  in  one  of  several  ways.  One  of  these  is  through  an 
assignment  statement  such  as  X = 10;  Another  is  through 
GET  statements  such  as  GET  LIST  (X);  or  GET  COPY  DATA. 
If  the  variables  are  a function  of  temperature  then  they 
should  be  initialized  by  calling  subroutine  READRK  as  fol- 
lows: 


CALL  READRK  (I,  IND,  Dl,  D2,  D3.  D4,  NI, 
INFILE); 


where: 


I is  an  index  (FIXED  BINARY). 

IND  the  name  of  the  table  contain- 

ing values  for  the  independent 
temperature  variable  (FLOAT 
BINARY  array  of  length  NI+1) 
<°R). 

Dl,  D2,  D3,  D4  are  variable  names  of  tables 

for  the  first,  second,  third, 
and  fourth  temperature  depen- 
dent variables  (FLOAT  BI- 
NARY numbers  of  length  NI), 

NI  is  a number  identifying  the 

length  of  the  tables  (FIXED 
BINARY),  and 

INFILE  is  the  data  file  from  which  the 

tables  are  to  be  read  (this 
variable  is  declared  as  a 
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FILE  in  the  declare  state- 
ments above.  If  the  input  is 
to  be  read  from  data  cards 
then  INFILE  is  SYSIN,  (FILE)). 

Data  cards  required  for  READRK  are  as  follows: 

Two  comment  cards,  followed  by  NI  cards,  each 
containing  the  appropriate  values  for: 

IND  in  columns  1 to  14 

D1  in  columns  15  to  29 


D2  in  columns  30  to  44 


All  in  F-Format 


D3  in  columns  45  to  59 


D4  in  columns  60  to  75 


Thermal  properties  of  each  material  must  be  initial- 
ized by  a call  to  READRK,  where: 

I is  a unique  number  from  1 

to  9 identifying  the  material, 

IND  is  temperature  (°R), 

D1  is  the  product  of  density  and 

specific  heat  (Btu/ft3.°R), 

D2  is  thermal  conductivity 

(Btu/lt  • h * °R), 


D3  is  surface  emissivity,  and 

D4  is  surface  absorptivity. 

Thermal  contact  conductances  must  also  be  input  via 
READRK,  where: 
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IND 

D1 

D2 

D3  and  D4 


is  a number  from  1 to  9 
identifying  the  conductance 
table  (this  number  must  be 
different  from  that  used  for 
material  property  tables 
and  other  conductance  tables), 

is  temperature  (°R), 

is  a dummy  variable  or  an- 
other variable  with  depen- 
dence on  temperature, 

is  contact  conductance 
(Btu/ft2  *h‘°R), 

are  dummy  or  other  varia- 
bles. 


Output  from  READRK  includes  the  two  comment  cards 
and  numerical  values  assigned  to  each  argument. 

5.  CALL  RSADTM  (I,  IND,  Dl,  D2,  D3,  D4,  D5, 

NL  INFILE);  reads  time  dependent  variables,  where: 

I is  an  index  defining  the 

READTM  statement  (FIXED 
BINARY), 

IND  is  an  array  containing  values 

for  the  independent  time  varia- 
ble (FLOAT  BINARY  array  of 
length  NI-<  1)  (seconds), 

Dl,  D2,  D3,  D4,  are  arrays  of  the  first,  sec- 

DS  ond,  third,  fourth,  and  fifth 

dependent  variables  (FLOAT 
BINARY  arrays  of  length  NI), 


u 
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is  the  length  of  the  above 
arrays  (FIXED  BINARY), 
and 
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INFILE  is  the  data  file  from  which 

above  arrays  are  to  be  read 
(this  variable  is  declared  as 
a FILE  in  the  declare  state- 
ments. If  input  is  read  from 
data  cards  then  INFILE  is 
SYSIN  (FILE)). 

Data  cards  required  for  READTM  are  as  follows: 

Two  comment  cards,  followed  by  NI  cards,  each  con- 
taining values  for: 


IND  in  columns 

1 to  9 

D1  in  columns 

10  to  19 

D2  in  columns 

20  to  29 

D3  in  columns 

30  to  39 

D4  in  columns 

40  to  49 

D5  in  columns 

50  to  59 

All  in  F-Format 


6.  CALL  READAB  — Subroutine  READAB  per- 
forms two  basic  functions.  First  it  computes  the  enthalpy  of 
the  ablative  material  in  its  solid  state.  This  is  accomplished 
by  trapezoidal  integration  of  the  capacitance /density  versus 
temperature  curve.  The  second  operation  of  READAB  is  to 
read  the  thermochemistry  data  supplied  by  the  EST  program 
(Ref,  2).  The  subroutine  reads  this  thermochemistry  data, 
changes  it  to  the  appropriate  units,  and  adapts  it  to  the  form 
required  by  subroutine'  ABAR03.  For  a Lewis  number  not 
equal  to  1 the  conversion  process  is  as  follows.  The  EST 

program  provides  values  of  gas  temperature  (T  ),  gas 

rg,  w 
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enthalpy  at  the  surface  (Hrg>  w),  the  enthalpy  of  the  gas  con* 
stitutes  at  the__boundary  layer  edge  evaluated  at  the  wall 

the  chemical  activity  function 

all  as  tabular  functions  of  gas 

pressure  (P^)  and  the  ablation  to  diffusion  mass  loss  ratio 
(0).  The  EST  data  is  expressed  in  calories,  grams,  and 
iC.  READAB  has  as  a parameter  the  conversion  factor  re- 
quired to  covert  these  values  to  Btu,  Ibm,  and_°R.  Further- 
more, the  routine  computes  log^  0,  Hrgw,  w,  and 
"QCHEM"  as  functions  of  and  Tw,  where: 


temperature  (H^  w),  and 


N 

"QCHEM"  = E <Z. 

i=l  l,e 


) AH. 


j.  A t u 


Notice  that  qchem  = Pu  <Cm)b  ’’QCHEM".  The  values  of  Tw 
used  in  this  table  are  specified  in  the  parameters  of  READAB. 

For  a Lewis  number  equal  to  1,  READAB  converts 
ESI'  values  as  follows.  The  EST  program  provides  values  of 
Tw  and  Hrgt  w as  functions  of  Pr  and  0.  The  data  is  con- 
verted to  units  of  Btu.  Ibm,  and  ^R.  In  addition  READAB 
also  computes  log10  0.  Hrg(W,  and  "QCHEM"  as  functions  of 
PL  and  Tw,  where  values  oi  Tw  are  specified  by  parameters 
READAB  and  "QCHEM"  = 0 <Ha>  w - Hrg  w).  Notice  that  in 
this  case  "QCHEM"  is  a misnomer  since  it  is  really  a por- 
tion of  qconv  and  (Eq.  (8)),  and  cannot  be  related  to 

%hemaltme- 


In  its  output  READAB  prints  the  ablative  material 
enthalpy  (Ha>w)  Eq.  (9),  as  a function  of  temperature  (Tw). 

It  also  prints  both  the  dimensionally  converted  input  data 
from  EST  as  functions  of  and  0,  and  the  reformated  data 
as  functions  of  P.  and  Tw.  Occasionally  one  may  not  wish 
to  have  these  tables  printed  in  the  output  as  they  can  be  quite 
lengthy.  If  so,  the  following  statement  is  placed  in  the  main 
program: 


DCL  AIRPRNT  INIT(O)  FIXED  BIN  EXT; 
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Subroutine  EEADAB  is  invoked  by  the  ceiling  se- 
quence: 

CALL  READAB  (#PL»  #BETA,  LEW#, 
RHO,  RHOCP,  TEMP  TWMIN,  TWMAX, 
TWINC,  CONV.  IN  FILE); 

where: 


#PL  specifies  the  number  of  pres- 

sure entries  being  supplied  by 
the  EST  program  (FIXED  BIN*, 

#BETA  specifies  the  number  of  ablation 

to  diffusion  mass  loss  ratio  03) 
values  being  supplied  by  the  EST 
program  at  each  value  of  pres- 
sure (FIXED  BINARY^ 

LEW#  is  the  Lewis  number  (FLOAT 

BINARY), 

RHO  is  the  surface  material  density 

(FLOAT  BINARY)  (lbm/ft3)t 

RHOCP,  TEMP  is  a list  of  density  times  specific 

heat  values  for  the  surface  mate- 
rial (from  which  values  of  HfttW 
are  calculated)  (the  values  cor- 
respond to  temperatures  in  the 
list  TEMP  (FLOAT  BINARY) 
(Btu/ftS,  °R);  the  length  of 
array  TEMP  is  one  greater 
than  array  RHOCP), 

TWMIN,  TWMAX,  are  constants  that  specify  the  be- 

and  TWINC  ginning  value  of  temperature,  the 

ending  value  of  temperature,  and 
the  increment  of  temperature  be- 
tween these  limits  for  the  variable 


! ! 
j ; 
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Tw,  which  it  used  in  convert- 
ing EST  data  from  fun  tion  of 
PL  and  0 to  functions  of  PL  and 


CONV  is  the  conversion  factor  re- 

quired to  convert  input  data  to 
°R.  If  data  is  in  6 K and  cal/g, 
then  CONV  =1.8;  otherwise 
CONV  =1.0  and  input  data 
have  units  of  °R  and  Btu/lbm 
(FLOAT  BINARY),  and 

INFILE  is  the  data  file  from  which  input 

to  READAB  is  to  be  read  (if  in- 
put is  read  from  data  cards  then 
INFILE  is  SYSIN  (FILE). 


READAB  automatically  assumes  an  enthalpy  reference  tem- 
perature (T ref)  of  536®  R.  If  a different  value  is  desired  the 
user  may  declare  TREF  as  FLOAT  BIN  EXT  in  the  main  pro- 
gram and  assign  it  the  appropriate  value  in  °R. 


In  addition  to  the  above  arguments  READAB  also  re- 
ceives information  from  data  cards  or  a data  tape.  These 
data  consist  of:  two  commentary  cards  followed  by  chemi- 
cal data  tables  obtained  from  the  EST  program.  The  data 
are  formatted  as  follows:  First  there  are  #BETA  cards  de- 
fining chemical  data  at  a given  local  pressure  (P.  ) as  a 
function  of  &.  This  set  of  cards  is  followed  by  subsequent 
sets  defining  chemical  data  at  other  local  pressure  levels. 
Hence  for  each  local  pressure  there  are  #BETA  cards  and 
there  are  #PL  pressures,  i.  e. , #PL*#BETA  data  cards  in 
the  table.  #BETA  and  #PL  must  be  greater  than  1. 


The  data  included  on  each  card  depend  on  the  Levis 
number.  If  LEW#  = 1 then  each  card  contains  the  following, 
in  F-format: 
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Units 


I 
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I 

I 

8 

D 
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0 

0 

0 

0 

0 

0 

0 

0 
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PL 

1 to  13 

atmospheres 

ft 

14  to  26 

none 

Tw 

27  to  39 

°C  or  °R 

Hrg,  w 

40  to  52 

cal/gm  or  Btu/lbm 

For  nonunity  Lewis  numbers  each  card  contains  the 
following,  in  F -format: 


Column 

Units 

PL 

1 to  13 

atmospheres 

ft 

14  to  26 

non* 

Tw 

27  to  39 

°C  or  °R 

Hrg,  w 

40  to  52 

cal/g  or  Btu/lbm 

Se.w 

53  to  65 

cal/g  or  Btu/lbm 

N 

T <Z.  - Z.  ) AH. 

i.e  i.w  i,  w 

66  to  78 

cal/g  or  Btu/lbm 

It  is  noted  here  that  the  above  described  variables 
should  not  be  of  direct  concern  to  the  user.  They  do  net 
enter  into  the  calling  sequence  of  any  of  the  subroutines  and 
their  allocation  and  subsequent  use  of  controlled  internally 
by  READAB  and  ABAR03.  Their  description  is  included 
only  for  completeness. 

Both  subroutines  (ABAR03  and  READAB)  make  use 
of  a number  EXTERNAL  variables  that  have  not  yet  been 
described.  These  variables  are  essentially  the  storage 
locations  for  the  thermochemical  data  (supplied  by  the  EST 
program)  as  needed  by  ABAR03  to  perform  the  ablation 
analysis.  These  variables  are; 
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PL 


BETA 


TW 


TWI 


QCHEM 


is  a one -dimensional  array  (or  list) 
of  pressure  values  (Pl,)  given  by  the 
EST  program  in  atmospheres,  and  is 
one  of  the  independent  variables  for 
the  property  data  (during  actual  use, 
this  variable  is  converted  to  the  base 
10  logarithm  of  the  pressure  value), 

is  a list  of  ablation  to  diffusion  mass 
loss  ratio  values  (9)  supplied  by  the 
EST  program  and  is  the  other  inde- 
pendent variable  for  the  given  prop- 
erty values  (during  actual  use  this 
variable  contains  the  natural  loga- 
rithm of  j8). 

is  a two-dimensional  array  of  values 
specifying  absolute  temperature  (°R) 
(each  row  of  this  dependent  variable 
is  a list  of  values  for  which  the  pres- 
sure is  constant  and  given  by  the  cor- 
responding value  in  PL), 

is  a one-dimensional  array  of  tem- 
peratures (°R)  specified  by  TWMIN, 
TWMAX,  and  TWINC  as  described 
above  (this  list  is  the  independent 
variable  used  with  PL  to  describe 
QCHEM,  HW,  and  HEW. 


is  a two-dimensional  array  of  either 
" N 


E (Z,  - Z.  ) AH. 

i=1  i.e  i,w  i. 


w 


+ ^<Ha.w-Hrg.w> 

for  a Lewis  number  not  equal  to  1, 
or 


^ ^a,  w ” **rg,  w) 
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for  a Lewis  number  equal  to  1,  con- 
verted by  READAB  to  depend  on  PL 
a ad  TWl  (Btu/lb), 

HW  is  a two-dimensional  array  of  wall 

enthalpy  (Hrg,  w)  converted  by 
READAB  to  uepend  on  PL  and  TWl 
(Btu/lb),  and 

HEW  is  a table  of  enthalpy  of  gases  at  the 

boundary  layer  edge  (H^  w)  evaluated 
at  Tw,  converted  to  depend  on  PL  and 
TWl  by  READAB  (Btu/lb). 

To  summarise  the  function  of  READAB  with  respect 
to  the  thermochemical  data,  the  following  sequence  is  listed; 


^ 1 . JThe  EST  program  is  run  to  produce  output  of 

Hrg.  w»  Tw.  He(W  and  the  chemical  activity  function: 


N 

z 

i=l 


<Zi.e-Zi.w)4Hi.w 


as  bivariant  tables  depending  on  local  pressure  and  abla- 
tion to  diffusion  mass  loss  ratio  0. 

2.  The  above  EST  output  is  rearranged  so  that  the 
pressure  and  0 values  are  ascending  in  magnitude  and  are 
then  reformatted  for  input  to  READAB  according  to  the  fol- 
lowing format: 


Format 

Columns 

Pl 

F-format 

1 to  13 

0 

F-  format 

14  to  26 

T 
* w 

F-format 

27  to  39 

Hrg,  w 

F-format 

40  to  52 
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H*tW  F- format  53  to  65 

N 

t izi  . “ z«  J AH,  «,  F-format  86  to  78 

M®  *|W  1|W  \ 

(There  la  available  a simple  reformatting  utility  program 
which  accepts  EST  format  and  produces  punched  output  in 
the  above  described  APL  format.  Note  that  for  a Lewis 
number  equal  to  1 the  last  two  items  may  be  omitted. ) 

3,  The  reformatted  EST  output  is  placed  in  the  in 
put  stream  of  the  CLAD  job  and  is  read  by  subroutine 
READAB. 


4.  READAB  applies  the  conversion  factor  CONV 
and  prints  out  the  converted  input.  This  printed  output  will 
contain,  in  addition  to  the  data  supplied,  the  calculated 
values  of  "QCHEM". 

5.  The  variables  Hr«  w and  "QCHEM"  are  then 
converted  to  depend  on  PL  and  Tw.  The  results  of  this  con- 
version are  then  printed  out. 

6.  Subroutine  ABAR03  uses  the  converted  data, 
stored  in  the  appropriate  EXTERNAL  variables  by  READAB, 
to  solve  the  required  surface  energy  balances. 

Since  the  thermochemical  data  are  supplied  as  tables, 
extensive  use  of  interpolation  is  required  to  obtain  the  pre- 
cise values  needed  in  ABAR03.  This  interpolation,  while 
not  of  direct  concern  to  the  user,  is  done  in  a logarithmic 
sense  with  regard  to  pressure  (P^)  and  (6).  When  these 
variables  are  required  for  use  in  interpolation,  their  loga- 
rithms are  first  taken  and  then  the  linear  interpolation  is 
performed.  The  antilog  of  this  interpolated  result  is  then 
supplied  as  the  required  value.  This  completed  the  discus- 
sion of  READAB. 
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7.  CALL  TRREED  and  CALL  REDCP3  - Geomet- 
ric data  necessary  to  describe  each  capacitor  in  the  thermal 
network  are  input  to  the  CLAD  program  in  one  of  two  ways. 
For  evenly  spaced,  rectangular  elements  the  user  can  call 
subroutine  THREED  to  automatically  divide  the  model  into 
elements  and  made  the  appropriate  thermal  connections. 

For  nonrectangular  elements  and  for  unevenly  spaced  rec- 
tangular elements  the  user  can  use  subroutine  REDCP3  to 
read  the  geometry  data.  These  two  subroutines  are  used 
as  follows: 

a.  CALL  THREED  - Subroutine  THREED 
divides  the  analytical  model  into  rectangular  elements  simi- 
lar to  Fig.  3a,  The  body  is  divided  into  finite  volume  sub- 
surface elements  and  zero  volume  elements  on  the  irradiated 
surface.  All  nodes  are  located  at  the  center  of  each  element. 
The  length  of  the  test  section  is  specified  for  the  x,  y,  and  z 
directions  with  each  length  divided  into  a specified  number 
of  elements.  In  the  z direction  the  model  may  be  made  up 
of  several  different  materials  separated  by  contact  conduc- 
tances. Each  material  and  contact  conductance  is  assigned 
an  identification  number  corresponding  to  the  number  (I)  in 
READRK.  The  largest  permissible  identification  number  is 
9.  However,  each  material  and  contact  conductance  may  be 
used  more  than  once,  which  permits  an  unlimited  number  of 
layers  in  the  Z direction, 

THREED  also  contains  an  option  w'hmh  allows 
for  a slightly  different  element  break-up  as  shown  in  Fig. 

3b.  In  this  configuration  elements  adjacent  to  the  xz  and  yz 
planes  are  half  as  large  as  the  interior  elements  and  have 
their  nodes  located  on  the  planes.  This  configuration  may 
be  particularly  useful  when  analyzing  stationary  laser  beams 
or  other  arrangements  having  two  planes  of  symmetry.  In 
both  models  elements  are  numbered  sequentially  starting  at 
x = 0.  y = 0,  z = 0 and  going  first  in  the  z direction,  then  x, 
then  y (as  shown  in  Fig.  3b).  In  addition  to  dividing  the 
model  into  elements,  THREED  also  computes  the  volume  of 
each  element  and  the  area  to  length  ratio  between  adjacent 
elements.  These  values  are  stored  in  the  array,  XCAP,  and 
passed  on  to  subroutine  CAPCN3  (see  Number  13  following). 
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CALL  THREED  (J,  LX,  NX,  LY,  NY,  NMZ, 

LZ,  NZ,  IPROP,  ICONT,  FEL,  SURF); 

where: 

J is  an  index  identifying  the 

THREED  statement  (FIXED  BI- 
. NARY),  '• 

LX,  NX,  LY,  NY  is  the  model  length  and  number 

of  nodes  in  the  x and  y direc- 
• tions  (lengths  LX  and  LY  are 

float  binary  numbers  in  feet, 
while  NX  and  NY  are  fixed  deci- 
mal values  not  greater  than  10), 

NMZ  is  the  number  of  material  layers 

in  the  z direction  (a  fixed  decimal 
number  from  1 to  “), 

LZ  is  an  array  defining  the  length  of 

each  material  layer  in  the  z di- 
rection starting  with  the  irra- 
diated material  (LZ  is  a float  bi- 
nary array  of  size  NMZ  with 
units  of  feet), 

NZ  is  an  array  cf  length  NMZ  defin- 

ing the  number  of  elements  in 
each  material  layer  (FIXED 
DECIMAL), 

IPROP  is  an  array  of  length  NMZ  identi- 

fying the  material  property  tables 
for  each  layer  (FIXED  DECIMAL), 
where  IPROP(l)  is  the  property  in- 
dex of  the  material  at  the  irra- 
diated surface  and  IPROP (2)  is  the 
property  index  of  the  adjacent  mate 
rial,  and  so  forth  (values  assigned 
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ICONT 


FEL 


to  IPROP  must  correspond  to  the 
appropriate  I in  READRK), 

is  an  array  of  length  (NMZ-1) 
identifying  tables  defining  con- 
tact conductivity  between  mate- 
rial layers,  starting  with  ICONT (1) 
for  the  conductance  between  the 
irradiated  material  (IPROP(l))  and 
the  adjacent  layer  (lPROP(2)) 
(FIXED  DECIMAL), 

is  the  index  number  of  the  ele- 
ment located  at  (0,  0,  0)  in  Fig.  3 
(FIXED  DECIMAL),  and 


SURF  is  a flag  (if  SURF  = 0. 0 the  model 

is  as  shown  in  Fig.  3a;  however, 
if  SURF  = 1.0  then  the  model  is 
as  shown  in  Fig.  3b  (FLOAT  BI- 
NARY). 

The  output  from  THREED  is  in  a tabular  format  with 
the  following  headings: 


II  VOLUME  MATERIAL  #CON  12  A/L  13  or  A A/L  14  A/L  T5  A/L 


These  headings  are  followed  by  a table  of  correspond- 
ing values.  The  headings  stand  for: 


VOLUME 


MATERIAL 


is  the  element  in  question, 

is  the  volume  of  element  II  (ft3), 

is  a single  digit  number  identifying 
the  material  of  element  II  (if  II  has 
a contact  conductance  connecting  it 
to  12  then  MATERIAL  is  a three 
digit  number,  the  first  digit 
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identifying  material  (II),  the  sec* 
ond  digit  identifying  the  contact  con- 
ductivity, and  the  third  digit  identi- 
fying material  (12)), 

#CON  specifies  the  number  of  conduction 

paths  from  II  to  elements  with  iden- 
tifiers greater  than  II, 

* 

12  is  the  element  to  which  II  conducts, 

A/L  is  the  area  to  length  ratio  of  element 

12  (feet)  (if  there  is  a contact  resis- 
tance between  II  and  12  then  A/L  is  the 
the  area  to  length  ratio  between  II 
and  the  contact  interface  (feet), 

13  or  A if  II  .d  12  have  a contact  resistance 

between  them,  then  13  is  the  area  of 
contact  interface  (ft2),  otherwise  13 
is  the  second  element  to  which  II  con- 
ducts, 

A/L  if  II  and  12  have  a contact  resistance 

between  them  then  A/L  is  the  area  to 
length  ratio  between  the  interface  and 
12  (feet),  otherwise  A/L  is  the  area 
to  length  ratio  between  elements  II 
and  13  (feet), 

14  is  the  element  to  which  II  conducts, 

A/L  is  the  area  to  length  ratio  between  II 

and  14  (feet), 

15  is  the  element  to  which  II  conducts, 
and 

A/L  is  the  area  to  length  ratio  between  II 

and  15  (feet). 
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The  reader  should  note  that,  by  convention, 
the  element  located  below  II  must  be  12.  However  if  there 
is  no  element  below  II  then  12  is  a side  element.  Also  note 
that  the  printout  contains  only  connections  to  elements 
located  below,  to  the  right,  and  behind  element  II. 

b.  CALL  REDCP3  (INFILE)  - Subroutine 
REDCP3  read§  the  geometric  data  necessary  to  describe 
each  capacitor  in  the  thermal  network.  The  only  argument 
to  REDCP3  is  INFILE  which  is  the  name  of  the  file  contain- 
ing the  data.  The  data  list  consists  of:  two  comment  cards 
followed  by  two  data  cards  for  each  capacitor  in  the  network 
to  be  described.  The  data  defining  each  capacitor  is  the 
same  as  discussed  in  the  printout  of  subroutine  THREED. 

The  necessary  card  format  and  units  are  listed 

below: 


Data 

Column 

Format 

Units 

Capacitor  number  (ID 

0 to  4 

Fixed 

Volume 

5 to  19 

E -format 

CO 

First  data 

Material 

20  to  29 

F-format 

card 

#CON 

30  to  39 

F -format 

12 

40  to  49 

F-format 

A/L 

50  to  59 

E -format 

feet 

13 

60  to  69 

F format 

A/L 

70  to  79 

E -format 

feet 

14 

0 to  9 

F-format 

A/L 

10  to  19 

E -format 

feet 

Second  data 

15 

20  to  29 

F-format 

card 

A/L 

30  to  39 

E -format 

feet 

- 48  - 


0 

0 


APPLIED  PHYSICS  LABORATORY 

•nvw  Mmm.  nA*iyj+m 


S 


I 

I 

I 


By  convention,  12  is  always  the  element  adja- 
cent *o  II  in  the  direction  of  ablation  (if  there  is  such  an 
element).  In  addition,  elements  in  an  ablating  string  must 
be  numbered  consecutively  starting  with  the  element  at  the 
ablating  surface. 

8.  CALL  SET  (START,  STOP.  TEMPIN,  L, 
INDEX,  TEMP);  — Subroutine  SET  specifies  the  start  and 
end  times  of  the  computer  analysis  and  sets  all  elements  to 
their  initial  temperature.  The  arguments  of  SET  are: 


START,  STOP 


TEMPIN 


start  computation  at  time  START 
and  stop  at  STOP  (FLOAT  BINARY) 
(seconds), 

all  elements  are  initialized  to  tem- 
perature TEMPIN  (FLOAT  BINARY) 
(°R). 

the  number  of  elements  to  be  initial- 
ized to  temperatures  other  than 
TEMPIN  (FIXED  BINARY),  and 

arrays  of  length  L defining  the  in- 
dex and  temperature  of  elements  not 
initialized  to  TEMPIN  (FIXED  BI- 
NARY and  FLOAT  BINARY),  where 
TEMP  has  units  of  °R. 


9.  II:  IMAX  = DECIDE  (TIMS,  IOX,  TIMED;  - 
This  statement  defines  the  maximum  intensity  (IMAX)  of  the 
Gaussian  beam  at  a specific  time,  TIMS: 


INDEX,  TEMP 


IMAX 


TIMS 


is  the  maximum  beam  intensity  at 
time,  TIMS  (FLOAT  BINARY  EX- 
TERNAL) (Btu/ft2  -s), 

is  the  present  time  (FLOAT  BINARY 
EXTERNAL)  (seconds). 
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is  an  array  of  maximum  intensity 
as  a function  of  time  (FLOAT  BI- 
NARY) (Btu/ft^  • s),  and 


TIME! 


is  the  corresponding  time  where: 
TIMEI(l)  is  the  number  of  entries 
in  the  table.  TIMEK2)  is  the  time 
corresponding  to  IOX(l).  etc. 
(FLOAT  BINARY)  (seconds). 


If  the  surface  is  not  irradiated  by  a Gaussian  beam 
then  the  above  statement  may  be  replaced  by  a series  of 
statements  defining  the  heating  rate  to  each  surface  node 
as  a function  of  time.  For  example: 


II:  Q(J) 


where: 


Q(J) 


DECIDE  (TIMS.  QJ,  TIMEJ); 


is  the  heating  rate  to  element  J 
(FLOAT  BINARY  EXTERNAL  CON- 
TROLLED) (Btu/h), 


, . | 
(.  * 


TIMS 


TIMEJ 


is  the  present  time  (FLOAT  BINARY 
EXTERNAL)  (seconds). 

is  an  array  defining  heating  rate  to 
element  J as  a function  of  time 
(FLOAT  BINARY)  (Btu/h),  and 

is  the  corresponding  time  (FLOAT 
BINARY')  (seconds)  where  TIMEJ(l) 

= length  of  table  and  the  remaining 
values  define  times  for  values  given 
in  QJ. 


10.  CALL  GAUSBM  (XO,  TIMEG,  XINIT,  RBEAM, 
XN,  TIMEP);  — GAUSBM  is  called  if  the  surface  is  being 
irradiated  by  a Gaussian  beam.  The  beam  may  be  stationary 
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or  may  move  along  the  x axis  as  a (Unction  of  time.  In 
either  case  the  beam  center  is  assumed  to  be  on  the  x axis. 
Within  GAUSBM  the  heating  rate  to  each  surface  element  is 
computed  from  the  equation: 

Q(J)  = IMAX  * AREA  * exp(-XN  * XY2/RBEAM2)  . (19) 

where  JMAX  is  the  maximum  beam  intensity  at  time  (TIMS), 
AREA  is  the  surface  area  of  element  J,  XY  is  the  distance 
from  the  beam  center  to  the  element  J,  XN  is  a constant, 
and  RBEAM  is  the  beam  radius. 

The  arguments  to  GAUSBM  are: 

XO  is  an  array  defining  the  distance 

from  x = 0 to  the  beam  center  mea- 
sured along  the  x axis  as  a function 
of  time  (FLOAT  BINARY)  (feet), 

TIMEG  is  the  corresponding  array  of  time 

(FLOAT  BINARY)  (seconds) 
TIMEG(l)  = length  of  array  XO, 

XINIT  is  the  initial  location  of  the  beam 

center  along  the  x axis  measured 
from  x = 0 (FLOAT  BINARY)  (feet), 

RBEAM  is  the  radius  at  which  the  intensity 

is  exp-XN  * imAX  (FLOAT  BINARY) 
(feet), 

XN  is  the  constant  in  Eq.  (19)  (usually 

either  1.0  or  2.0)  (FLOAT  BINARY), 
and 


TIMEP 


is  an  array  listing  desired  printout 
times  of  surface  intensity,  usually 
the  same  as  specified  in  WRITE 
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At  times  specified  by  T1MEP  subroutine  GAUSBM 
prints  the  following: 


BEAM  DEFINITION 
CENTERED  ATX*  FT  IMAX  * 


BTU/S  FT  **  2 


IRRADIATED  ELEMENT  AREA  = FT  **  2 

INTENSITY  AT  EACH  ELEMENT  (BTU/S  * FT  **  2) 


I ( ) 


I ( 


) = 


I ( 


) 


etc. . etc. 


If  the  elements  are  not  divided  such  that  each  ele- 
ment has  the  same  surface  area,  then  the  IRRADIATED 
ELEMENT  AREA  given  above  is  the  area  of  elements  not 
adjacent  to  either  the  x or  y axes. 

11.  CALL  QFUSON;  - QFUSON  computes  the  latent 
heat  of  fusion  remaining  as  an  element  is  heated  or  cooled 

at  its  melting  point.  Values  of  latent  heat  of  fusion  (QMELT). 
melt  temperature  (TMELT)  and  density  (RO)  for  each  mate- 
rial are  obtained  from  the  main  program  as  EXTERNAL 
variables. 

12.  CALL  ABAH03  (I,  #SURF,  CYL,  LAMBDA. 

TSP,  CM_CH,  LEW#,  VIEW,  FACTOR,  TLIM,  QLIM, 

ITLIM,  TIMEP,  IPLOT);  — ABAR03  computes  the  tempera- 
ture of  each  zero  volume  surface  node  using  the  energy  bal- 
ance given  in  Eq.  (1).  The  equation  is  solved  iteratively  us- 
ing temperature  dependent  values  of  surface  absorptivity, 
material  thermal  conductivity,  emissivity,  and  chemical 
energy.  ABAR03  also  removes  material  from  ablating  and 
vaporizing  elements,  and  adjusts  conduction  paths  to  adja- 
cent elements. 

Arguments  of  ABAR03  are: 


I 


is  an  index  identifying  the  ABAR03  state- 
ment (FIXED  BINARY), 


is  the  number  of  surface  elements  being 
aerodynamical ly  heated  in  the  call  state- 
ment (FIXED  BINARY), 

is  a flag  for  cylindrical  or  rectangular 
geometry  (FIXED  BINARY), 

Input:  0 for  rectangular  elements 

I for  heating  on  the  outer  surface 
of  cylindrical  elements 
-1  for  heating  on  the  inner  surface 
of  cylindrical  elements, 

is  the  "blowing  parameter"  in  Eq.  (13) 
(FLOAT  BINARY)  (a  complete  discussion 
of  LAMBDA  is  presented  in  Ref,  3 — the 
"classical"  value  is  LAMBDA  = 0.  5;  how- 
ever, Ref,  3 recommends  that  LAMBDA  = 

0.  5 for  laminar  flow  and  LAMBDA  - 0.  4 
for  turbulent  data), 

is  the  radiation  space  temperature  (FLOAT 
BINARY)  (°R), 

is  the  ratio  of  mass  transfer  Stanton  number 
to  heat  transfer  Stanton  number  (FLOAT 
BINARY)  (if  CM__CH  is  negative  then  the 
mechanical  erosion  is  included  as  presented 
in  Eq.  (18),  and  the  Putts  and  Bartlett  blow- 
ing equations  are  used  - Eqs.  (14)  through 
(17)), 


is  the  Lewis  number  (FLOAT  BINARY), 


is  the  view  factor  of  incoming  laser  irra- 
diation (FLOAT  BINARY)  (note  this  value 
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applies  only  to  laser  irradiation  and  not  to 
qrad;  see  Eq.  (D), 

is  a lumping  factor  (if  the  volume  of  chemi- 
cally ablating  or  vaporising  element  drops 
below  this  fraction  of  its  original  volume, 
then  it  is  lumped  with  the  element  below 
{FLOAT  BINARY)). 

is  the  temperature  convergence  limit  of 
energy  balance  Eq.  (1)  (FLOAT  BINARY) 
<°F). 

is  the  heat  convergence  limit  of  energy 
balance  (FLOAT  BINARY)  (Btu/ft2-s), 

is  the  maximum  number  of  iterations  al- 
lowed for  convergence  of  energy  balance 
to  within  the  two  tolerance  values  listed 
above  (FIXED  BINARY);  if  this  limit  is 
exceeded  then  the  following  message  is 
printed  and  the  program  continues: 


ITERATION  ON  SURFACE  NODE  # HAS  EXCEEDED TSURF  = 

AND  DELTA  TSURF  - . HEAT  IMBALANCE  » , 

TIMEP  is  an  array  of  printout  data  (FIXED  BINARY) 

(seconds). 

I PLOT  is  a flag  for  plotting  data  (FIXED  BINARY) 

If  0 then  plots  are  not  desired  (if  nonzero 
then  many  of  the  values  computed  in  ABAR03 
are  stored  in  the  temperature  vector  and 
may  be  plotted  by  calling  subroutine  PLOTTER; 


i 
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for  more  information  on  the  use  of  thie  rou- 
tine the  reader  la  referred  to  Ref.  IS).  The 
following  values  are  stored  in  the  tempera- 
ture array  T: 


T(IPLOT) 


TOPLOT+l)  * total  mass  loss  (lbm/ft2), 
TUPLOT+2;  - mass  loss  rate  (lbm/ft2  *s), 


TftPLOT  3)  35  accumulated  error  in  the  surface  energy 

balance  (Btu/ft2), 


T(IPLOT+4) 

TUPLOT+5) 

T(IPLOT+6) 

T(IPLOT+7) 


* current  error  in  the  surface  energy  balance 
<Btu/ft2  • s), 

= current  recession  rate  (in/s). 

* total  recession  (inches). 

= convective  heating  rate  (PeUe  (CH)b 

(Hr  - He  w))  when  the  Lewis  number  does 
not  equal'  i;  and  when  the  Lewis  number 
equals  1 then  TUPLOT+7)  is  he  firjst  term 
in  Eq.  (8),  i.e.  <PeUe  (CH>b  (Hr  - H » 
(Btu/ft2  *s). 


T(IPLOT+8)  = accumulated  conduction  away  from  the  sur- 
face (Btu/ft2), 

T(IPLOT+9)  = current  heating  rate  from  the  "QCHEM  ma 
term  (Btu/ft2  »s)  where:  for  a Lewis  number 
not  equal  to  1. 


Ref.  13.  R.  K,  Frazer,  Subroutine  for  Obtaining 
CalComp  Plots  of  SHTP  Variables,  " APL/JHU  EM-4496, 
13  December  1972. 
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TdPLOT+lO) 

TdPLOT+11) 

TOPLOTY  12) 
TdPLOT+13) 
TdPLOT+14) 


T(IPLOT+15) 

TdPLOT+16) 

TdPLOT+17) 

TdPLOTU8) 


and  for  a Lewis  number  equal  to  1, 

■^CHEM".4(Ha(W-Hrgw). 

> accumulated  heating  from  the  '^QCHEM"  * mt 
term  (Btu/ft2), 

* current  radiative  energy  flux  absorbed  at 
the  surface,  including  laser  radiation  and 
radiation  from  a radiating  body  (Btu/ft2  • a), 

* current  heat  radiating  from  the  surface 
(Btu/ft2*s), 

* current  conduction  loss  rate  from  the  sur- 
face <Btu/«2.s), 

* current  rate  of  sensible  energy  required  to 
heat  the  ablating  material  from  the  node 
temperature  (Tg  in  Fig.  1)  to  the  surface 
temperature  (Tw)  (Btu/ft2*s), 

* surface  node  number, 

* current  value  of  surface  node  area  (ft2). 

* current  radius  of  surface  node  (feet),  and 

e accumulated  sensible  energy  required  to 
heat  ablating  material  from  T,  to  Tw 
(Btu/ftS).  2 


Since  the  above  information  is  stored  in  temperature 
array  (T)  it  is  necessary  that  IPLOT  be  greater  than  the 
largest  element  number  in  the  analysis.  Also,  be  sure  to 
set  LASCAP  to  a large  enough  value  in  statement  3. 
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In  addition  to  its  arguments.  ABAR03  also  receives 
information  in  the  form  of  data  cards.  Data  input  in  this 
manner  are  the  surface  area  and  radius  of  each  surface  ele- 
ment aerodynamically  heated  in  this  call  statement.  There- 
fore. there  will  be  #SURF  cards,  each  containing:  the  sur- 
face element  number,  the  surface  area  (ft2),  and  the  surface 
radius  (feet).  If  the  elements  are  not  cylindrical  a radius  of 
0. 0 may  be  input.  There  is  no  column  or  format  restriction 
on  these  input  data. 

13.  CALL  CAPCN3;  - CAPCN3  computes  capaci- 
tance of  each  element  and  heat  conduction  between  elements 
using  geometry  values  placed  in  array  XCAP  by  THREED  or 
REDCP3,  No  arguments  are  required  for  CAPCN3,  since 
all  values  are  obtained  from  previously  computed  EXTER- 
NAL variables. 

14.  CALL  WRITE  (TIMEP);  - Subroutine  WRITE 
prints  out  element  temperatures  at  times  specified  in  array 
TIMEP. 

TIMEP  is  an  array  of  desired  printout 

times  (FLOAT  BINARY)  (seconds) 
(if  TIMEP(l)  is  a negative,  non- 
*ero  number,  the  output  will  be 
printed  for  each  time  step). 

15.  CALL  STEP  (MINSTP,  MAXSTP):  - Subroutine 
STEP  computes  the  stability  time  step  and  the  element  tem- 
peratures at  the  next  step  in  time, 

MINSTP  is  the  minimum  allowable  time 

step  (FLOAT  BINARY)  (seconds), 
and 

MAXSTP  is  the  maximum  allowable  time 

step  (FLOAT  BINARY)  (seconds). 

If  the  time  step  is  greater  than  MAXSTP  then  it  will  be  set 
to  MAXSTP  and  the  program  continues.  However,  if  the 
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time  step  is  less  than  MINSTP,  the  following  message  is 
printed  and  the  program  is  terminated. 


TIME  STEP  LESS  THAN  MINIMUM  ALLOWED 

CRITICAL  INDEX  NO.  DFLTIM  = 

If  subroutine  WRITE  has  printed  during  the  present  time 
step  then  STEP  prints  the  following: 


CRITICAL  INDEX  NO.  = TIME  STEP  = SEC. 

$ & $ # $ j{<  # # >jt  >!<  >!< )(« sfc  >js  >{<  >;<>(<  >j<  >!<  & 


The  row  of  asterisks  signifies  the  end  of  printout  for  the 
present  time  step.  When  the  end  time  has  been  reached 
(i.  e.  time  STOP  as  defined  in  the  argument  of  SET),  the 
following  message  is  printed  and  the  program  is  terminated: 

THE  TIME  IS SECONDS,  AND  THE  NUMBER  OF 

STEPS  EXECUTED  WAS . 

16,  GO  TO  II;  — Sends  control  of  the  program  back 
to  II,  which  is  the  first  statement  in  the  time  loop  (see  Step 
9).  Note  that  the  loop  is  exited  via  the  routine  step  and  is 
not  infinite. 


17.  END  NAME;  — This  statement  signifies  the 
end  of  MAIN  program  NAME. 
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4.  SAMPLE  PROBLEM 


Each  phase  of  the  CLAD  program  has  been  checked 
against  available  analytical  and/or  experimental  data.  The 
main  program  used  to  check  ablation  computations  in  CLAD 
is  presented  here  as  a sample  problem.  A schematic  of 
the  analytical  model  is  shown  in  Fig.  4.  The  model  con- 
sists of  a cylindrical,  ATJ-S  blast  tube,  convectively  heated 
on  the  inner  surface.  From  Ref.  14  the  enthalpy  based  con- 
vective heat  transfer  coefficient  is  3.  96  lb/ft2  • s and  the  re- 
covery temperature  is  4880°R.  The  chemical  ablation  data 
for  ATJ-S  graphite  was  obtained  from  the  EST  computer 
program  (Ref.  3)  using  a reference  temperature  of  536°R. 
Hence,  the  recovery  enthalpy  above  536°R  is  the  enthalpy  of 
air  at  4880°  R minus  the  enthalpy  of  air  at  536,  which  is 
1416  - 129  or  1287  Btu/lb. 

The  cylinder  is  1.  4 inches  thick  and  has  an  inner 
radius  of  0,435  inch.  The  analysis  begins  at  0.2  second 
with  the  model  temperature  distribution  shown  in  Fig.  4. 

The  main  program  (ABMAIN)  used  in  the  analysis  is 
presented  in  Appendix  B.  The  reader  should  refer  to  Appen- 
dix B as  the  statements  of  ABMAIN  are  discussed  below. 
ABMAIN  is  divided  into  sections  headed  by  comment  cards: 

/*  DECLARE  ENTRY  STATEMENTS  */ 

Under  this  heading  are  entry  points  for  all  subrou- 
tines called  in  ABMAIN. 

/*  DECLARE  VARIABLES  */ 


These  statements  declare  variable  names  used  in 
call  statement  arguments.  These  variables  are  discussed 
later,  as  they  are  encountered  in  the  CALL  statements. 


Temperature  and  Ablation  Histories  for  a Proposed  Blast 
Tube  at  ARC,  " APL/JHU  EM-4547,  7 January  1974. 
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/*  ALLOCATE  AND  INITIALIZE  VARIABLES  */ 


The  allocate  statement  assigns  space  to  controlled 
variables  RO,  QVAPOR.  and  TVAP,  while  the  GET  COPY 
DATA  statement  assigns  initial  temperatures  to  elements  2 
through  11.  Subroutine  STORE  assigns  an  array  size  of 
LASCAP  to  allocatable  variables  used  in  the  laser  heating 
subroutines.  The  remaining  statements  assign  printout 
times  to  variable  TIMEP. 


/*  READ  INPUT  DATA  */ 


In  this  section  values  are  assigned  to  variables 
which  have  not  been  previously  initialized  in  the  above  de- 
clare statements.  The  first  few  statements  read  and  print 
values  for  heat  of  vaporization  (QVAPOR),  vapor  tempera- 
ture (TVAP),  and  material  density  (RO)  of  the  first  mate- 
rial. These  statements  are  followed  by  call  statements 
which  read  input  values. 


READRK  assigns  thermal  properties  to  material  1 
as  a function  of  temperature. 


READTM  assigns  time  dependent  values  to  recovery 
enthalpy  (HR1),  radiation  heating  rate  (QRAD),  cold  wall  en- 
thalpy based  heat  transfer  coefficient  (RHOCH1)  and  local 
pressure  (PLTM1). 


Since  the  Lewis  number  equals  1,  READAB  reads 
local  pressure  and  P dependent  values  of  _wall  tempera- 
ture (Tw)  and  reacted  gas  wall  enthalpy  (Hrg>w).  Using 
these  values  READAB _computes  "QCHElVr  i.e., 

Ij8  (Ha>w  - Hrg>w)l  and  H w as  functions  of  log1Q  (PL) 
mperature  (TWL 


^ ®a,  w ~ **rg,  w^ J anc* 
and  wall  temperature  (Tw 


REDCP3  reads  the  volume  of  each  element  and  the 
conduction  paths  between  elements. 
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/*  INITIALIZE  MODEL  TEMPERATURE  */ 


CALL  SET  specifies  the  analysis  to  run  from  0. 2 
(TZERO)  to  0.  5 (ENDTIME)  second.  SET  also  initializes 
elements  12  through  18  to  540°R  (TINIT)  and  assigns  ele- 
ments 2 through  11  the  initial  temperatures  read  in  GET 
COPY  DATA. 


/*  STATEMENTS  IN  THE  TIME  STEP  LOOP  */ 

Subroutine  ABAR03  computes  model  surface  tem- 
peratures. ablation  rates,  heat  of  vaporization,  and  chemi- 
cal heat  of  reaction,  and  adjusts  conduction  paths  to  ablat- 
ing elements.  The  first  argument  indicates  this  is  the  first 
call  to  ABAR03.  There  is  one  (#SURF)  surface  node  and  it 
is  located  on  the  inside  of  a cylinder  (CYL  = -1).  The  blow- 
ing parameter  constant  is  0.4  (LAMBDA)  and  the  surface 
radiates  to  a 0.0°  R space  temperature  (TSP).  The  ratio  of 
mass  transfer  to  heat  transfer  is  1,  as  is  the  Lewis  number. 
If  there  were  laser  irradiation,  it  would  impinge  on  the  sur- 
face with  a view  factor  of  1 (VIEW).  When  half  (FACTOR) 
an  element  is  ablated  it  is  then  lumped  with  the  adjacent  ele- 
ment. Convergence  criteria  for  the  surface  temperature 
iteration  are  ± 0. 1°  F (TLIM)  and  ± 0.  5 Btu/ft2  (QLIM).  The 
iteration  stops  after  10  (ITLIM)  iterations.  Printout  times 
are  specified  in  TIMEP  and  there  are  no  plots  (PLTCODE  = OX 

CALL  CAPCN3  computes  heating  rates  between  ele- 
ments and  has  no  arguments. 

CALL  WRITE  specifies  printout  times  (TIMEP). 

CALL  STEP  indicates  the  critical  time  step  is  not  to 
be  less  than  0.  005  second  (MINSTP)  or  greater  than  0.  5 sec- 
ond (MAXSTP). 

TO  TO  IMX:  sends  control  of  the  program  back  to 
the  beginning  of  the  time  loop. 

END  ABMAIN:  signifies  the  end  of  MAIN  procedure. 
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The  main  program  is  followed  by  data  cards  which 
are  listed  in  Appendix  B in  the  order  they  are  required  in 
the  MAIN  procedure.  Data  cards  define  the  following  varia- 
bles: 


1.  The  number  of  elements  to  be  assigned  initial 
temperatures  and  a table  of  temperatures  and 
corresponding  element  numbers. 

2.  Heat  of  vaporization,  vapor  temperature,  and 
density  of  material  number  1.  Zero  vapor  tem- 
perature indicates  the  material  is  ablating  and 
not  vaporizing. 

3.  Temperature  dependent  values  of  capacitance, 
conductivity,  emissivity,  and  absorptivity  for 
material  number  1. 

4.  Time  dependent  values  of  recovery  enthalpy,  in- 
put radiation,  enthalpy  based  heat  transfer  coef- 
ficient, and  local  pressure  at  the  surface. 

5.  Pressure  and  B dependent  values  of  wall  tem- 
perature (Tw)  and  gas  enthalpy  at  the  wall 
(Hrg,  w)  f°r  the  ablative  material  in  air. 

6.  Geometry  values  defining  each  element,  its  vol- 
ume, its  material,  the  number  of  elements  con- 
nected to  it  via  conduction,  and  the  adjacent 
conduction  elements  and  corresponding  area  to 
length  ratio. 

7.  The  surface  area  and  radius  of  element  1. 

This  ends  the  sample  problem  computer  program. 
The  reader  should  realize  this  program  has  been  written  in 
a general  form  which  allows  its  application  to  a great  many 
problems  simply  by  changing  a few  data  cards.  This  makes 
parametric  studies  very  easy  and  eliminates  recompilation 
of  the  MAIN  program  for  each  run. 
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Printout  resulting  from  the  above  program  is  shown 
in  Appendix  C.  The  first  six  and  one  half  pages  consist  of 
input  data  for  the  GET  COPY  DATA.  GET  LIST,  READRK, 
READTM,  READAB,  REDCPS,  and  ABAR 03 statements  of 
the  MAIN  program  (see  above  discussion  of  input  data). 

This  is  followed  by  a series  of  values  computed  for  each 
time  step.  At  a typical  time  step  (0.231?  second)  the  first 
call  to  ABAR03  prints  the  following  values:  recovery  en- 
thalpy = 1287  Btu /lb,  cold -wall  enthalpy  based  neat  transfer 
coefficient  = 3.960  lb/ft2  • s,  local  pressure  = 100  atmo- 
spheres. ratio  of  mass  transfer  to  heat  transfer  = 1,  Lewis 
number  = 1,  blowing  coefficient  - 0.4,  and  radiation  space 
temperature  = 0.  0001°R.  Node  number  1 has  a "convective 
heating  rate"  of  2925.  3 Btu/ft2  *s  and  a cumulative  "convec- 
tive heating  rate"  of  90.  2 Btu/ft2.  When  the  Lewis  number 
is  equal  to  1,  as  in  this  case,  the  "convective  heating  rate" 
is  a misnomer  and  is  really  equal  to  the  first Jterm  in  the 
right  side  of  Eq.  (8),  i.e. . Pe Ue  (CH)b  (Hr  - Hrg>w).  Con- 
tinuing with  the  printout  in  Appendix  C,  there  is  no  incoming 
radiation  and  the  mass  loss  rate  is  0.  6484  lb/ft2  • s with  a 
total  mass  loss  of  0.  0210  lb/ft2.  The  accumulated  error  in 
the  surface  energy  balance  is  -0.00164  Btu/ft2.  Node  num- 
ber 1 is  at  4381.  8°F,  The  heating  rate  into  the  surface  due 
to  "QCHEM"  is  740  Btu/ft2  • s and  the  total  heating  due  to 
"CHEM"  is  24.4  Btu/ft2.  There  is  no  radiation  from  the 
surface.  The  recession  rate  is  0.0701  in/s  with  a total  re- 
cession of  0.00227  inch.  The  heat  imbalance  for  the  pres- 
ent time  step  is  0. 142  Btu/ft2  *s.  The  mass  loss  coefficient 
(BETA)  equals  0. 175.  Heat  conduction  away  from  the  sur- 
face is  3218.  7 Btu/ft2  • s with  an  accumulated  heat  conduc- 
tion loss  of  104.6  Btu/ft2.  The  sensible  energy  to  increase 
the  surface  material  from  T2  to  Ti  is  446.  2 Btu/ft2  • s with 
an  accumulated  value  of  4.  91  Btu/ft2.  The  ratio  of  heat 
transfer  coefficient  with  blowing  to  no  blowing  is  0.  9359. 

The  wall  gases  have  an  enthalpy  (Hrg|W)  of  497.6  Btu/lb 
and  the  mass  flow  Stanton  number  is  3.  706.  The  surface 
energy  balance  data  is  followed  by  temperature  data.  The 
temperature  of  each  element  is  printed  in  degrees  Fahren- 
heit. The  final  printout  for  this  time  step  is  the  critical  in- 
dex number  (6)  and  the  critical  time  step  (0.0112  second). 
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Surface  temperature  histories  from  this  analysis 
are  presented  in  Fig.  5.  Notice  the  sawtooth  effect  caused 
by  nodal  lumping  as  the  surface  ablates.  This  sawtooth  ef- 
fect can  be  greatly  reduced  by  using  more  elements  near 
the  surface.  With  smaller  elements  the  true  surface  tem- 
perature will  fall  between  the  peaks  and  valleys  of  the  CLAD 
curve  presented  in  Fig.  5.  Also  shown  in  Fig.  5 are  results 
from  the  one-dimensional  CMA  program  (Ref.  3).  CMA  re- 
sults are  slightly  lower  because  it  uses  a finite  volume  sur- 
face node  while  CLAD  uses  a zero  volume  node.  Also  note 
that  because  CMA  drops  nodes  from  the  back  surface  it  does 
not  have  any  sawtoothing  effect.  The  apparently  better  back 
node  dropping  technique  is  very  difficult  to  implement  in  a 
three-dimensional  node-dropping  analysis  and  consequently 
was  not  used  in  CLAD. 
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5.  CONCLUSIONS  AND  RECOMMENDATIONS 


A comprehensive  computer  program  has  been  written 
to  evaluate  material  damage  to  laser  irradiated  flight  bodies. 
The  CLAD  program  includes  many  features  such  as  aerody- 
namic heating,  radiation  relief,  and  temperature  dependent 
thermal  properties  that  are  also  contained  in  similar  laser 
da  mage  computer  programs.  However,  CLAD  contains  sev- 
eral unique  features  which  are  not  generally  included  in  the 
othe*  programs  (e.  g. , three-dimensional  material  removal, 
chemical  reactions  and  vaporization  at  the  surface,  and  in- 
depth  melting).  The  combination  of  all  these  features  into  a 
single  program  provides  an  analytical  tool  that  should  be 
very  useful  in  analyzing  laser  heated  flight  vehicles. 

There  are  also  several  areas  where  CLAD  can  be 
expanded  and  improved.  For  example  it  could  be  expanded 
to  include  melt  removal  by  aerodynamic  shear  forces  at  the 
surface  or  to  include  in-depth  chemical  reactions  which 
occur  in  charring  ablators. 


[ 


The  reader  should  also  realize  that  the  CLAD  pro- 
gram was  written  for  continuous  wave  laser  applications, 
which  have  relatively  low  power  densities  compared  with 
pulsed  lasers.  At  the  high  intensities  of  a pulsed  laser, 
ordinary  thermal  analysis  no  longer  holds.  Consequently 
the  CLAD  program  will  have  to  be  extensively  modified  for 
application  with  pulsed  lasers.  These  modifications  are 
planned  for  the  future. 

Even  without  the  above  modifications,  CLAD  is  still 
a powerful  analytical  tool  for  a wide  range  of  laser  heated 
material  analyses. 
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DERIVATION  OF  EQUATION  USED  TO  COMPUTE 
CHEMICAL  HEATING  RATE 


The  components  of  the  chemical  heating  term  in  Eq. 
(1)  are  shown  schematically  in  Fig.  A-l.  In  Fig.  A-l. 
gases  diffuse  through  the  boundary  layer  at  a rate  of  (md) 
and  have  an  average  total  enthalpy  H^g>w.  At  the  surface 
some  of  the  gases  react  with  the  ablating  material  <ma)  and 
release  chemical  energy  qcbem.  The  resulting  gases  then 
mix  with  the  unreacted  diffusion  gases  and  form  a gaseous 
mixture  whose  average  total  enthalpy  is  Hrg  w.  The  re- 
sulting mixture  leaves  the  surface  via  diffusion  and  convec- 
tion at  a rate  of  (mrg).  The  rate  of  heat  released  at  the 
surface  is: 


q . = m . H , m n 

Mchem  d dg,  w a a,w 


+ m H 


m H . 
rg  rg,  w 


(A-l) 


Notice  that  all  enthalpies  in  Eq.  (A-l)  are  total  enthalpies 
(i.  e. , they  are  the  sum  of  the  sensible  and  chemical  en- 
thalpies). 


Obtaining  Eq.  (A-l)  in  the  form  shown  in  Eq.  (5)  of 
the  text  is  accomplished  as  follows.  Since  there  is  no  mass 
build-up  at  the  surface  and  the  mass  flow  out  (mrg)  equals 
the  mass  flow  in  (ma  + m^).  Substituting  for  mrg*in  Eq. 
(A-l)  and  rearranging  qchem  results  in: 


q . = m , (H  . - H ) + m 

^chem  d dg,  w rg,w  i 


(H  - H ).  (A-2) 
a,  w rg,  w 


However,  ma  can  be  defined  by  a parameter,  Q, 
where  ma  B d PeUe  <Cm)b.  For  an  informative  discussion 
of  the  physical  meaning  of  d.  see  Ref.  14. 
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la  Eq.  (A-2)  the  total  diffusioa  mass  flow  rate  (m^) 
arriving  at  the  wall  is  the  sum  of  the  flow  rates  of  individ- 
ual species,  each  arriving  at  enthaly  H|  w and,  after  chemi- 
cal reaction,  leaving  at  Hrg  w.  Hence,  the  energy  term 

N 


(H, 


^dg,w  " Hrg,w)  in  Eq*  <A_1)  becomes  S mi(H. 


-H  ). 
w rcr.  w 


Furthermore,  according  to  Ref.  5,  diffusion  mass  transfer 
can  be  related  to  concentration  coefficients  at  the  boundary 
layer  edge  (Zi>e)  and  at  the  wall  (Z^  w)  as  follows: 


N 

z 

i=l 


m. 

1 


p U 
e e 


<C  L 
m b 


N 

Z 

i=l 


(Z.  - Z.  ) . 

1,  e i,  w 


(A-3) 


In  addition,  the  enthalpy  difference  (H^  w - Hrg  w) 
is  simply  the  chemical  enthalpy  change  of  each  specie 

^Hi,  w^* 

Hence,  Eq.  (A-2)  becomes: 


=p  U 
e e 


(C  >K 
m b 


N 

z 

i=l 


(Z.  -Z.  ) AH.  +/3  (H  -H  ),  (A-4) 

i.e  i,w  l,  w a,  w rg,  w 


which  is  in  the  form  presented  in  Eq.  (5). 


- 73  - 


APPLIED  PHYSICS  LABORATORY 

•dm*  |N(M.  MAAtlAMO 


0 

0 

0 

u 

y 

y 

o 

y 


y 

i 

!' 

i - 


Appendix  B 

MAIN  PROGRAM  FOR  SAMPLE  PROBLEM 

This  appendix  presents  the  complete  program  used 
in  analyzing  the  sample  problem  discussed  in  Section  4. 
The  program  contains  the  following  comment  statements 
which  describe  its  overall  organization. 

/*  DECLARE  ENTRY  STATEMENTS  */ 

/*  DECLARE  VARIABLES  */ 

/*  ALLOCATE  AND  INITIALIZE  VARIABLES  */ 

/*  READ  INPUT  DATA  */ 

/*  INITIALIZE  MODEL  TEMPERATURE  */ 

/*  STATEMENTS  IN  THE  TIME  STEP  LOOP  */ 

ABARO 

CAPCN3 

WRITE 

STEP 

The  program  is  followed  by  input  data. 
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Appendix  C 

SAMPLE  PROBLEM  PRINTOUT 


A portion  of  the  sample  problem  printout  is  pre- 
sented in  this  appendix.  The  first  page  contains  data  read 
into  the  main  procedure  via  GET  statements  and  READ  sub- 
routines. The  next  five  pages  contain  thermochemical 
tables  for  ATJ-S  graphite  and  air,  and  are  followed  by  two 
tables  defining  the  model  geometry.  The  remaining  pages 
present  typical  printouts  at  0.2  second.  0.221  second,  and 
0.  232  second.  The  time -dependent  output  is  divided  into 
aerodynamic  heating  data,  surface  energy  balance  informa- 
tion, thermal  capacitor  temperatures,  and  critical  time 
step  data. 
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NOMENCLATURE 


area 

absorptivity  of  the  surface  evaluated  at  Tw 

heat  transfer  Stanton  number  (qConv^e^e  AH) 

mass  transfer  Stanton  number  (m^/PgUg  AZ) 

specific  heat 

natural  logarithm  base 

radiation  view  factor 

the  sensible  enthalpies  above  Tref  of  an 
ablative  material  evaluated  at  Tw  and 
T2.  respectively 

sensible  recovery  enthalpy 

sensible  gas  enthalpy  at  the  wall 

chemical  enthalpy  change  of  each  diffusing 
specie  i due  to  reactions  with  the  abla- 
tive material.  AH^  w is_evaluated  at 
Tw,  and  equals  Hi/W  - Hrg>w 

total  enthalpy  above  Trej  (i.  e. , H = sensi- 
ble + chemical  enthalpies) 

total  enthalpy  above  Tref  of  solid  ablative 
materials  evaluated  at  Tw 

average  total  enthalpy  above  Tref  of  gases 
diffusing  to  the  surface,  evaluated  at 

Tw 

average  total  enthalpy  above  Tref  of  gases 
at  the  thermal  temperature  boundary 
layer  edge,  evaluated  at  Tw 

total  recovery  enthalpy  above  Tre£ 
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Hj  w total  enthalpy  above  Tref  of  gaseous 

specie  i diffusing  to  the  surface, 
evaluated  at  Tw 

w average  total  enthalpy  (above  Tref)  of  the 

reacted  gases  at  the  surface,  evaluated 
at  Tw 

Ki_^  thermal  conductivity  between  elements  i 

and  j 

Lj_j  distance  from  element  i to  element  j 

ma  ablative  mass  loss  rate  per  unit  area 

m^  mass  flow  rate  per  unit  area  of  all  gases 

diffusing  to  the  surface 

ni|  mass  flow  rate  per  unit  area  of  gaseous 

species  i diffusing  to  the  surface 

mmech  mechanical  erosion  mass  loss 


m 


rg 


^chem 


^cond 

^conv 

^d 


mass  flow  rate  per  unit  area  of  all  reacted 
gases  leaving  the  surface  via  diffusion 
and  mass  transfer 

number  of  species  diffusing  to  the  wall 

local  pressure 

chemical  heating  rate  per  unit  area  due  to 
ablative  material  reacting  at  the  sur- 
face 

chemical  heating  rate  per  unit  area 
^chem  ^d  + 3a) 

conduction  heat  transfer  rate  per  unit  area 

convection  heating  rate  per  unit  area 

chemical  heating  rate  per  unit  area  due  to 
diffusing  gases  reacting  at  the  surface 

total  heating  rate  per  unit  area  due  to 
laser  heating  and/or  other  radiative 
heating 
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qlas 
^rad  in 

^rad  out 

qsens 
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vap 

ref 


incident  laser  heating  rate  per  unit  area 

heating  rate  per  unit  area  due  to  radiation 
to  the  surface 

heat  loss  rate  per  unit  area  due  to  radia- 
tion from  the  surface 

heating  rate  per  unit  area  required  to 
raise  the  sensible  energy  of  the  solid 
ablative  material  from  T2  to  Tw 

heat  of  vaporization  (Btu/lb) 

reference  tempei  ature  at  which  chemical 
and  sensible  enthalpies  are  zero  (Tref 
is  often  set  at  536°R) 

space  temperature,  to  which  the  surface 
radiates 


Tw,  Tg.  Tg  temperature  at  the  wall  and  at  interior 

nodes  2 and  3 


U, 


velocity  parallel  to  the  surface  at  the 
boundary  layer  edge 


V 

VIEW 
zi,  e*  zi,w 


a 

0 
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element  volume 

view  factor  for  the  incident  laser  beam 


diffusion  mass  transfer  potentials  of 
specie  i evaluated  at  the  boundary- 
layer  edge  and  at  the  wall  respectively 
(these  potentials  are  discussed  in  Ref.  5) 

2Xm 

blowing  coefficient,  — - — — — 

Pe  e CH 


ablation  to  diffusion  mass  loss  ratio  de- 

m 


fined  as  # = 


P V (C  ). 
e e mb 


thermal  emissivity  of  the  surface  evaluated 
at  T, 
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Subscripts 

a 

b 

d 

dg 
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i 

rg 
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an  empirical  constant  known  as  the  "blow- 
ing parameter"  (usually  X = 0. 5 for 
laminar  flow  and  X = 0. 4 for  turbulent 
flow) 

material  density 

density  at  the  boundary  layer  edge 

Stefan-Boltsman  constant 

times  at  three  consecutive  time  steps 
(a,  b,  and  c) 


ablative  material 

corrected  for  blowing 

diffusion 

diffusion  gas 

boundary  layer  edge 

specie  i diffusing  to  the  wall 

radiating  gas 

recovery  condition 

wall 

from  node  2 to  node  5 
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CONTRACTORS 


Waahlngton,  D.  C. 
China  Lake,  Cal. 


San  Diago,  Cal. 
Crane,  ind. 


Washington,  D.  C. 
WPAFB,  Ohio 
WPAFB,  Ohio 
Edwards  AFB.  Cal. 

Arnold  AFS,  Tsnn. 
St.  Louia,  Mo, 

San  Antonio,  Texas 
Andrews  AFB.  Md. 

Maxwell  AFB,  Ala. 


Redstone  Arsenal,  Ala. 
Huntsville,  Ala. 

Aberdeen  Proving  Ground,  Md. 
Dover,  N.,1. 


Washington,  D.  C. 
Ft.  Meade,  Md. 


Washington,  D.  C. 
Washington,  D,  C. 

Hampton,  Va. 

Cleveland,  Ohio 
Edwards  AFB.  Cal. 
Morfett  Field,  Cal. 
College  Park,  Md. 

Greenbelt,  Md. 


Automation  Industries /Vitro  Labs.  Div.  Silver  Spring,  Md. 


Aarodyn.  Lab.  Llbr. 

406 

45 

Tech.  Library,  753 
Library 

Naval  Ammo.  Depot 


AFRDD 
F.  D.  Stull 
AVO-2 
RPMMM 
MKCO 

DSC/Reaearch.  AER 

PRA 

PDEG 

TDT 

SCTSP 

Tech.  Llbr. /DPSL 
Tech.  Library 


AMCPM-MDEM 
Chief,  Doc.  Section 
AMXBR-XM-SE 
STINFO 


CSR/ADD/Std.  Dist. 
C3/TDL 


RAPL  N.  Rekos  (Code  RLCJ 
M.  J,  Aueremanne 
(Code  SG) 

K.  F.  Rubert 
J.  Henry 

E.  Lexberg 

J.  Nugent  (Prop. ) 

F.  Pfyle  228.  ? 

J.  Waldo 


Library  - 252 


Report  Acquis.  1SA-U 


«li>we  attribution  of  the  dowtmmt  within  the  ApsMad  fbyUo  Laboratory  hai  bean  lead*  In  aceocdmct  with  a Wat  ee  Wa  Intha  APL  Tachntot  PubHwtlom  Qwup. 
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